#include "PT_table.h"

bool PT_Table::import_table(string filename)
{
    // 从文件中导入表
    fstream file;
    file.open(filename, ios::in);

    // 如果文件打开失败则返回 false
    if(!file.is_open())
    {
        cout << "Failed to open file: " << filename << endl;
        return false;
    }

    while(!file.eof())
    {
        string line;
        getline(file, line);

        // 若开头不是 PVTTABLE POINT 则跳过
        if(line.find("PVTTABLE POINT") == string::npos)
            continue;

        // 提取第一个括号之后的所有内容
        int start = line.find("(");
        string content = line.substr(start + 1);

        // 提取出所有使用逗号分隔的数字
        vector<double> numbers;
        int last = 0;
        for(int i = 0; i < content.size(); i++)
        {
            if(content[i] == ',' || content[i] == ')')
            {
                numbers.push_back(stod(content.substr(last, i - last)));
                last = i + 1;
            }
        }
        table_ori.push_back(numbers);

        // // 打印到控制台
        // for(auto num : numbers)
        //     cout << num << ", ";
        // cout << endl;
    }
    file.close();

    // 遍历所有 P T, 建立索引表
    for(auto& entry : table_ori)
    {
        double P = entry[0];
        double T = entry[1] + 273.15;

        // 若 P 不存在则插入
        P_index_map[P] = 0;
        // 若 T 不存在则插入
        T_index_map[T] = 0;
    }

    // 为每个 P T 建立索引
    int index = 0;
    for(auto& entry : P_index_map)
    {
        entry.second = index;
        index++;
    }
    index = 0;
    for(auto& entry : T_index_map)
    {
        entry.second = index;
        index++;
    }

    // // 打印信息到控制台
    // cout << "P_index_map: " << endl;
    // for(auto& entry : P_index_map)
    // {
    //     cout << entry.first << " -> " << entry.second << endl;
    // }
    // cout << "T_index_map: " << endl;
    // for(auto& entry : T_index_map)
    // {
    //     cout << entry.first << " -> " << entry.second << endl;
    // }

    // 初始化表格
    table.resize(P_index_map.size());
    for(auto& row : table)
    {
        row.resize(T_index_map.size());
    }

    // 将原始数据填入表格
    for(auto& entry : table_ori)
    {
        PT_Entry pt_entry;
        pt_entry.P = entry[0];
        pt_entry.T = entry[1] + 273.15;// 表格中的单位为摄氏度，转换为开尔文
        pt_entry.rho_gas = entry[2];
        pt_entry.rho_liq = entry[3];
        pt_entry.rs = entry[11];
        pt_entry.visc_gas = entry[13];
        pt_entry.visc_liq = entry[14];
        pt_entry.heat_gas = entry[16];
        pt_entry.heat_liq = entry[17];
        pt_entry.ent_gas = entry[19];
        pt_entry.ent_liq = entry[20];
        table[P_index_map[pt_entry.P]][T_index_map[pt_entry.T]] = pt_entry;
    }
    return true;
}

void PT_Table::save_to_mat(string filename)
{
    // 保存表格到文件
    fstream file;
    file.open(filename, ios::out);

    // 打印列名
    file << ",";
    for(auto& T_index : T_index_map)
    {
        file << T_index.first << ",";
    }
    file << endl;

    // 打印表格
    for(auto& P_index : P_index_map)
    {
        file << P_index.first << ",";
        for(auto& T_index : T_index_map)
        {
            auto& entry = table[P_index.second][T_index.second];
            file << entry.rho_gas << ",";
        }
        file << endl;
    }


    file.close();
}

// 使用双线性插值法，得到 P T 对应的密度和粘度
PT_Entry PT_Table::find(double P, double T)
{
    // 找到 P T 对应的四个点
    double P1 = 0, P2 = 0, T1 = 0, T2 = 0;

    // 使用 lower_bound 找到 P1 和 P2
    auto it = P_index_map.lower_bound(P);
    // 若没找到，说明 P 太大，直接使用最大值
    if(it == P_index_map.end())
    {
        P2 = (--it)->first;
        P1 = (--it)->first;
        P = P2;
    }
    // 若找到的是第一个，说明 P 太小，直接使用最小值
    else if(it == P_index_map.begin())
    {
        P1 = it->first;
        P2 = (++it)->first;
        P = P1;
    } 
    // 否则找到两个点
    else
    {
        P2 = it->first;
        P1 = (--it)->first;
    }

    // 使用 lower_bound 找到 T1 和 T2
    it = T_index_map.lower_bound(T);
    // 若没找到，说明 T 太大，直接使用最大值
    if(it == T_index_map.end())
    {
        T2 = (--it)->first;
        T1 = (--it)->first;
        T = T2;
    }
    // 若找到的是第一个，说明 T 太小，直接使用最小值
    else if(it == T_index_map.begin())
    {
        T1 = it->first;
        T2 = (++it)->first;
        T = T1;
    }
    // 否则找到两个点
    else
    {
        T2 = it->first;
        T1 = (--it)->first;
    }

    // 计算四个点的权重
    double w_P1_T1 = (P2 - P) * (T2 - T);
    double w_P2_T1 = (P - P1) * (T2 - T);
    double w_P1_T2 = (P2 - P) * (T - T1);
    double w_P2_T2 = (P - P1) * (T - T1);
    // 权重归一化
    double w_sum = w_P1_T1 + w_P2_T1 + w_P1_T2 + w_P2_T2;
    w_P1_T1 /= w_sum;
    w_P2_T1 /= w_sum;
    w_P1_T2 /= w_sum;
    w_P2_T2 /= w_sum;

    // 拿出四个点的密度和粘度
    PT_Entry pt_entry_P1_T1 = table[P_index_map[P1]][T_index_map[T1]];
    PT_Entry pt_entry_P2_T1 = table[P_index_map[P2]][T_index_map[T1]];
    PT_Entry pt_entry_P1_T2 = table[P_index_map[P1]][T_index_map[T2]];
    PT_Entry pt_entry_P2_T2 = table[P_index_map[P2]][T_index_map[T2]];

    // 计算插值结果
    PT_Entry result;
    result.P = P;
    result.T = T;
    result.rho_gas =    w_P1_T1 * pt_entry_P1_T1.rho_gas + 
                        w_P2_T1 * pt_entry_P2_T1.rho_gas + 
                        w_P1_T2 * pt_entry_P1_T2.rho_gas + 
                        w_P2_T2 * pt_entry_P2_T2.rho_gas;
    result.rho_liq =    w_P1_T1 * pt_entry_P1_T1.rho_liq +
                        w_P2_T1 * pt_entry_P2_T1.rho_liq +
                        w_P1_T2 * pt_entry_P1_T2.rho_liq +
                        w_P2_T2 * pt_entry_P2_T2.rho_liq;
    result.rs      =    w_P1_T1 * pt_entry_P1_T1.rs +
                        w_P2_T1 * pt_entry_P2_T1.rs +
                        w_P1_T2 * pt_entry_P1_T2.rs +
                        w_P2_T2 * pt_entry_P2_T2.rs;
    result.visc_gas =   w_P1_T1 * pt_entry_P1_T1.visc_gas +
                        w_P2_T1 * pt_entry_P2_T1.visc_gas +
                        w_P1_T2 * pt_entry_P1_T2.visc_gas +
                        w_P2_T2 * pt_entry_P2_T2.visc_gas;
    result.visc_liq =   w_P1_T1 * pt_entry_P1_T1.visc_liq +
                        w_P2_T1 * pt_entry_P2_T1.visc_liq +
                        w_P1_T2 * pt_entry_P1_T2.visc_liq +
                        w_P2_T2 * pt_entry_P2_T2.visc_liq;
    result.heat_gas =   w_P1_T1 * pt_entry_P1_T1.heat_gas +
                        w_P2_T1 * pt_entry_P2_T1.heat_gas +
                        w_P1_T2 * pt_entry_P1_T2.heat_gas +
                        w_P2_T2 * pt_entry_P2_T2.heat_gas;
    result.heat_liq =   w_P1_T1 * pt_entry_P1_T1.heat_liq +
                        w_P2_T1 * pt_entry_P2_T1.heat_liq +
                        w_P1_T2 * pt_entry_P1_T2.heat_liq +
                        w_P2_T2 * pt_entry_P2_T2.heat_liq;
    result.ent_gas =    w_P1_T1 * pt_entry_P1_T1.ent_gas +
                        w_P2_T1 * pt_entry_P2_T1.ent_gas +
                        w_P1_T2 * pt_entry_P1_T2.ent_gas +
                        w_P2_T2 * pt_entry_P2_T2.ent_gas;
    result.ent_liq =    w_P1_T1 * pt_entry_P1_T1.ent_liq +
                        w_P2_T1 * pt_entry_P2_T1.ent_liq +
                        w_P1_T2 * pt_entry_P1_T2.ent_liq +
                        w_P2_T2 * pt_entry_P2_T2.ent_liq;
    
    return result;
}